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Although anger is an important emotion that underlies much 
overt aggression at great social cost, little is known about how to 
quantify anger or to specify the relationship between anger and the 
overt behaviors that express it. This paper proposes a novel statistical 
model which provides both a metric for the intensity of anger and an 
approach to determining the quantitative relationship between anger 
intensity and the specific behaviors that it controls. From observed 
angry behaviors, we reconstruct the time course of the latent anger 
intensity and the linkage between anger intensity and the probabil- 
ity of each angry behavior. The data on which this analysis is based 
consist of observed tantrums had by 296 children in the Madison WI 
area during the period 1994-1996. For each tantrum, eight angry be- 
haviors were recorded as occurring or not within each consecutive 
30-second unit. So, the data can be characterized as a multivariate, 
binary, longitudinal (MBL) dataset with a latent variable (anger in- 
tensity) involved. Data such as these are common in biomedical, psy- 
chological and other areas of the medical and social sciences. Thus, 
the proposed modeling approach has broad applications. 

1. Introduction. Anger is an important emotion that can intrude into 
daily life as low intensity irritation. At higher intensity, it underlies overt 
aggression at great social cost. For example, it has long been established 
that anger can trigger partner abuse and assault [e.g., Burgess et al. (2001), 
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Jacobson et al. (1994), Schumacher et al. (2001)]. More recently, it has been 
found to play a role in incidents of road rage [e.g., Lupton (2002), Parker, 
Lajunen and Summala (2002)]. Anger at its most intense has been claimed 
to play a causal role in 35% of homicides [Curtis (1974)]. To approach this 
socially important but poorly understood phenomenon scientifically, we need 
to be able to measure and quantify anger intensity. Unfortunately, little is 
known about quantifying anger [cf., Fridja et al. (1992)]. The self-ratings 
of anger intensity on subjective rating scales (e.g., 1-10) used in classic 
emotion research are not open to verification. In contrast, people's behavior 
can reliably indicate the intensity of their anger. For instance, a grimace or 
a grunt suggests irritation, a shout indicates anger, and a screaming assault 
demonstrates rage. According to Sonnemans and Frijda (1994), the severity 
of angry action is one of the strongest predictors of the overall intensity of 
felt anger. Thus, an alternative model of anger would involve quantifying 
anger intensity based on a set of observable behaviors. 

In a typical episode, anger first rises and then falls. What little is known 
about this trajectory from self-report studies suggests that anger intensity 
rises rapidly and then declines slowly [e.g.. Beck and Fernandez (1998), 
Fridja et al. (1991), Tsytsarev and Grodnitzky (1995)]. Characterizing this 
trajectory is essential in understanding how different angry behaviors are 
distributed within an episode of anger. However, anger episodes vary widely 
in duration [see Potegal, Kosorok and Davidson (1996) for a review]. One 
resulting complication is that the trajectory of anger might vary system- 
atically with duration. Thus, a complete model must specify not only the 
relationship between underlying anger and the overt angry behaviors, but 
also how the trajectory of anger varies as a function of episode duration. 
These are the major goals of this paper. 

A problem in behavior-based modeling is that angry behaviors in adults 
are highly idiosyncratic and therefore difficult to compare or collapse across 
subjects. Furthermore, adults tend to mask their emotions in public situa- 
tions and their angry behaviors are correspondingly difficult to observe. In 
contrast, the high frequency of young children's tantrums indicates that they 
tend not to mask their emotions [Underwood, Coie and Herbsman (1992)]. 
Also, the angry behaviors appearing in tantrums are stereotyped and simi- 
lar across children. They are easily observed and identified by parents. For 
these reasons, our study focuses on young children's tantrums. 

Earlier work suggests that tantrum behaviors, like the angry acts of 
adults, are differentially associated with low and high intensities of anger. For 
instance, stamping is associated with lower anger intensity, and screaming 
and shouting reflect higher anger intensity [Potegal and Davidson (2003), 
Potegal, Kosorok and Davidson (2003)]. In this context, it is of consider- 
able interest that higher intensity anger behaviors are most likely to occur 



MODELING THE TIME COURSE OF ANGER 



3 



relatively early in tantrums, while lower anger behaviors are more evenly dis- 
tributed [Potegal, Kosorok and Davidson (2003)]. These observations could 
be confirmed if the probability of lower anger behaviors did not change much 
with anger intensity, but the probability of higher anger behaviors were to 
increase strongly with anger intensity. Our analysis will determine if this is 
the case. 

The data forming the basis of this analysis are derived from written 
parental narratives of tantrums of 296 children ranging from 18 to 60 months 
old, which were collected in the Madison WI area during the period 1994- 
1996 by Potegal and colleagues [Potegal and Davidson (2003), Potegal, 
Kosorok and Davidson (2003)]. Parents described each of their child's tantrums 
in as much detail as possible, including the durations of individual events or 
sets of events, within the tantrum. Detailed written instructions, checklists 
and examples guided parents in writing their narratives. Coders then con- 
verted the written narratives into time versus behavior matrices in which 
time was partitioned into consecutive 30-second units and eight different 
angry behaviors were scored as occurring or not within each unit. These 
behaviors are stiffen limbs/arch back, shout, scream, stamp, push/pull, hit, 
kick and throw. They are denoted as xi,X2, ■ ■ ■ ,xs, respectively, in this pa- 
per. A tantrum begins with the first occurrence of one of these behaviors and 
it is over when all these behaviors disappear. This tantrum dataset can be 
characterized as a multivariate, binary, longitudinal (MBL) data associated 
with a latent variable, anger intensity. The observed angry behavior variables 
xi,X2, ■ ■ ■ ,xs are assumed to be driven by the (unobservable) anger inten- 
sity. Several indicators of data reliability, reported in Potegal and Davidson 
(2003), suggest that this dataset is reasonably reliable. 

To model MBL data, generalized linear modeling [McCullagh and Nelder 
(1989)] and generalized estimating equations (GEEs) procedures [Diggle, 
Liang and Zeger (1994), Liang and Zeger (1986)] are natural options. Con- 
ventional GEE procedures would assume a generalized linear model for mean 
responses, with anger duration, time, etc. as predictors. Several authors, 
including Lin and Carroll (2001) and Severini and Staniswalis (1994), gen- 
eralized the GEE method by including a nonparametric component in the 
mean responses in addition to a parametric component. However, these ex- 
isting methods may not be appropriate for the current data because they do 
not allow us to study the relationship between the observed angry behav- 
iors and the unobservable anger intensity. Wu and Zhang (2002) suggested 
a nonparametric procedure for modeling longitudinal data. But that proce- 
dure is for cases with univariate continuous responses, and it does not allow 
any latent predictors. 

In this paper, a statistical modeling approach is suggested for describing 
the tantrum data. Our proposed model has three major components. First, 
we assume that the latent anger intensity, here called momentary anger 
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(MA), follows a flexible parametric form over the time course of a tantrum. 
Although certain existing psychological research suggests an asymmetrical 
trajectory of MA with a rapid rise and slow decline, as mentioned above, 
our model must be capable of representing a range of possible trajectories. 
These include a high initial value (i.e., the anger intensity rises so rapidly 
that it appears instantaneous under our conditions of observation) followed 
by a slow decline over the entire tantrum, a symmetrical rise and fall, or a 
gradual increase over the entire tantrum. After an extensive comparison of 
various candidate functions, including polynomial functions and the Gamma 
function, the two-parameter Beta function is chosen for modeling MA (see 
Section 2 for a detailed description). Second, potential dependence of MA 
on tantrum duration is handled by allowing each of the two parameters of 
the Beta function to be a polynomial function of duration. Third, the rela- 
tionship between MA and the likelihood of each of the eight angry behaviors 
is assumed to follow a generalized polynomial model. Model parameters are 
estimated by a proposed iterative algorithm, which is a modified version of 
the conventional GEE algorithm. In that algorithm, all model parameters 
are grouped into several blocks to speed up computation. 

MBL data with one or more latent variables involved are quite popular in 
biomedical, psychological and other areas of medical and social sciences. For 
instance, observable sleeping status (i.e., deep sleep, light sleep or awake) of 
animals is believed to be affected by the unobservable circadian rhythm [cf., 
Qiu (2002), Qiu et al. (1999)]. The proposed modeling approach has broad 
applications in these and many other studies. 

The remainder of this article is organized as follows: In Section 2 our 
proposed model for describing the tantrum data is described in detail. In 
Section 3 we apply this approach to the tantrum data and report some 
results. Several remarks conclude the article in Section 4. Identifiability of 
the proposed model and model estimation are discussed in the Appendices. 

2. Proposed model. The tantrum dataset studied in this paper consists 
of 296 tantrums of widely varying duration. The shortest duration is only 
0.5 minutes, while the longest tantrum lasts for 39.5 minutes. All 296 dura- 
tions are summarized and displayed in Figure 1. In this plot, the durations 
are classified into the following six categories: 0.5-2 minutes, 2.5-4 minutes, 
4.5-10.5 minutes, 11-20 minutes, 20.5-30 minutes, and 30.5-39.5 minutes. 
For convenience in plotting, these categories are slightly modified to the class 
intervals: (0.25,2.25], (2.25,4.25], (4.25,10.75], (10.75,20.25], (20.25,30.25], 
and (30.25,39.75]. The densities of these class intervals, which are defined 
by the relative frequencies of durations in the intervals divided by the cor- 
responding lengths of the intervals, are displayed in this plot as a density 
histogram. So, the area of each bar denotes the probability that a randomly 
chosen duration would be in the corresponding interval. In this plot, the 
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integer above each bar denotes the frequency of the durations that are in- 
cluded in the corresponding interval. It can be seen that only about 10% of 
the tantrums are longer than 11 minutes. 

One major goal of this research project is to describe temporal variation 
of MA and the relationship between MA and the eight angry behaviors in an 
entire tantrum episode. To this end, we first standardize observation times of 
angry behaviors along a to 1.0 time scale. That is, the standardized time, 
denoted as t, for an observation in a given tantrum is the actual time of 
that observation divided by the tantrum duration. In other words, the stan- 
dardized observation time t is the percentile of the duration of that tantrum 
episode at the given observation time. Let 7rk{t) denote the probability of 
the kth angry behavior at time t, for A; = 1, 2, . . . , 8, and MA(t) be the latent 
variable MA at t. Then, in this section, we discuss statistical modeling of 
MA{t) and the linkage between 7rk{t) and MA{t). 

As noted above, some existing psychological research suggests that MA{t) 
may increase rapidly at the beginning of the tantrum and then decline 
gradually. However, to avoid introducing a premature bias into determin- 
ing the trajectory of MA, candidate functions must be able to assume a 
range of possible shapes. As mentioned above, after searching quite exten- 
sively among many commonly used parametric families of curves defined on 
[0, 1] that could be used for describing MA(t), we found that the following 
two-parameter family of Beta curves can serve this purpose well: 

MA{t,r,s) = f~\l-ty^^ for t€ [0,1], 

where r,s > are two parameters. Figure 2 presents several Beta curves. It 
can be seen that the two parameters of the Beta function make it highly 
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Fig. 1. Density histogram of the durations of the tantrum dataset. 
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Fig. 2. Beta curves when parameters (r, s) take the values listed in the legend box. 



polymorphic. By changing values of these two parameters, it can assume 
many different shapes, including the three possible shapes of MA(t) men- 
tioned in Section 1. 

In the above expression of MA(t, r,s), the restriction that r and s are 
positive would be inconvenient for model estimation discussed in Appendix 
B. For this reason, we reparameterize that expression by replacing r by 
and s by e'', where a and b are two new parameters taking arbitrary values 
on the entire line R. After this reparameterization, MA(t,r, s) becomes 

(1) MA{t,a,b)=t^''-\l-tf-^ forte [0,1]. 

For simplicity of presentation, sometimes we write MA{t,a,b) as MA(t), 
which should not cause confusion. 

To allow possible dependence of MA(t,a, 6) on the tantrum duration, 
denoted as d, the two parameters a and b are assumed to be the following 
polynomials of d: 

(2) a = ao + aid + --- + amad^'', b = bo + bid + ■ ■ ■ + bn^^d""" , 

where and rrih are two non-negative integers. As shown in Figure 1, 
tantrum durations can change quite dramatically. In the psychological lit- 
erature, researchers doubt that the trajectory of MA(t) may depend on 
tantrum duration [cf., Fridja et al. (1992), Potegal, Kosorok and Davidson 
(1996)]. Models (1) and (2) have the flexibility to allow such dependence. 
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As mentioned in Section 1, the eight observable angry behavior variables 
xi,X2, ■ ■ ■ iXs are assumed to be driven by MA, each with its own linkage 
function. To model the linkage between xi,X2, ■ ■ ■ ,xq and MA, the following 
generalized linear model is used: for k = 1,2, ... ,8, 

(3) logit(7rfc(t, a, b, Cfc)) = cofc + cik MA{t, a, 6) + • • • + Cm,k MA™'^ {t, a, b), 

where logit(7rfc(t, a, b, Cfc)) = log(7rfc(t, a, b, Cfc)/(1 - 7rfc(t, a, b, Cfc))), a= (ao, 
ai,...,ar,ij',h = {bo,bi,...,bmJ,Ck = (cofc, cifc, . . . , 0^^^)', and are non- 
negative integers. In the above expression, we have made the coefficients a, b 
and Cfc explicit in the notation of Tr^it, a, b, c^), for convenience of discussion 
about model estimation. Model (3) allows the log-odds of the likelihood of 
each behavior variable to be a polynomial function of MA(t, a, 6). In Ap- 
pendix A, we show that parameters in models (l)-(3) are all identifiable. 

3. Analysis of the tantrum data. This section is organized in two parts. 
We first present a simulation study in Section 3.1 about the model estimation 
algorithm used in analyzing the tantrum data. Then, some results about the 
tantrum data are presented in Section 3.2. 

3.1. A simulation study. The proposed models (l)-(3) are estimated by 
a modified version of the generalized estimating equations (GEE) algorithm. 
In that modified version, instead of updating all parameters of the models 
c = (a', b', c']^ , . . . , Cg)' simultaneously in each iteration of the GEE algorithm, 
we suggest dividing c into g blocks, that is, c = (c'^, £'2, . . . , c^)', and then 
updating the parameters block-by-block in each iteration of the algorithm. 
By doing so, computation can be greatly reduced. Details of the GEE algo- 
rithm and our proposed version with blocked parameters are described in 
Appendix B. 

To check potential impact of the blocking scheme on parameter estima- 
tion by our modified version of the GEE algorithm [cf., expression (15) 
in Appendix B], we perform the following simulation study. It is assumed 
that there are = 200 subjects in the study. For the ith subject, 3 binary 
variables are observed at rii equally spaced time points tij £ [0, 1] , where 
j = 1,2, ... ,ni and i = 1, 2, . . . , A^. Therefore, durations of all tantrums are 
assumed to be 1 in this example. It is further assumed that = 5 when 
1 < i < 50, nj = 6 when 51 < i < 100, rij = 7 when 101 <i< 150, and = 8 
when 151 < i < 200. Probability of occurrence of the kth. binary variable at 
time tij is assumed to follow the model 

(4) logit(7rfc(ty)) = cok + cik(tij + l3tij) for /c = 1, 2, 3, 

where true values of parameters are set to be cqi = 0.5, C02 = 0.5, cqs = 
0,cii = C12 = C13 = 1, and (3 = —1. Obviously, model (4) is a special case of 
models (l)-(3) when = e'' = 2 in (1) and rrik = 1 in (3) for all k = 1, 2, 3. 
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Table 1 

This table presents the averaged estimates and the corresponding standard errors (in 
parentheses) of the parameters of model (4) by algorithm (15) in Appendix B with the 
three blocking schemes B-I, B-II and B-III. The results are based on 100 replicated 

simulations 



Parameters 


True values 


B-I 


B-II 


B-III 


coi 


0.5 


0.5273 (0.0185) 


0.5279 (0.0185) 


0.5279 (0.0185) 


Cll 


1 


0.8770 (0.0980) 


0.8694 (0.0981) 


0.8694 (0.0981) 


C02 


0.5 


0.5105 (0.0157) 


0.5132 (0.0159) 


0.5132 (0.0159) 


Cl2 


1 


0.9287 (0.0869) 


0.9168 (0.0878) 


0.9168 (0.0878) 


Cos 





-0.0074 (0.0156) 


-0.0101 (0.0153) 


-0.0101 (0.0153) 


CIS 


1 


1.1041 (0.0788) 


1.1152 (0.0773) 


1.1152 (0.0773) 


P 


-1 


-1.0170 (0.0224) 


-1.0061 (0.0211) 


-1.0061 (0.0211) 



In the modified GEE algorithm (15), we consider the following three block- 
ing schemes: 

B-I: g = l and ci = (cqi, cn, co2, ci2, cqs, C13, ^)', 
B-II: 5 = 2, ci = (coi,cii,co2,ci2,co3,ci3)', and C2 = (3, 
B-III: 5 = 4, ci = (coi,cii)', C2 = (co2,ci2)', C3 = (003,013)', and C4 = (3. 

Obviously, algorithm (15) with blocking scheme B-I is the same as the con- 
ventional GEE algorithm [cf., expression (14) in Appendix B]. Blocking 
scheme B-II divides all parameters into 2 blocks, while blocking scheme 
B-III divides the parameters into 4 blocks. 

In (15) initial values of (cq/cCi/j) are set to be their logistic regression 
estimates of model (4), for k = 1,2,3. Initial value of [3 is set to be the aver- 
age of the three logistic regression estimates of (3 obtained from model (4) 
when k = 1,2,3. From 100 replicated simulations, Table 1 presents the av- 
eraged estimates of all parameters by algorithm (15) and the corresponding 
standard errors (in parentheses), when the three blocking schemes defined 
above are used. From the table, we can see that (i) parameter estimates 
with different blocking schemes are almost identical, and (ii) algorithm (15) 
estimates the parameters reasonably well. 

3.2. Real data analysis. In this part we present some results about the 
tantrum data described in Section 2, using the modeling procedure (l)-(3). 
To estimate parameters in (2) and (3), the proposed modified version of 
the GEE algorithm (15) is used, after all parameters are grouped into the 
following 9 blocks: (a', b')', ci, . . . , cg. From Table 1, we know that different 
blocking schemes would have minimal effect on the performance of (15). 

To implement (15), we first need to choose a set of initial values c for the 
parameters. While there are no existing general guidelines for this purpose in 
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the GEE literature, we use the following procedure for choosing c^^^ First, 
in (2), we let a'^'^ = log(1.5), Sq'^^ = log(3), and the remaining components of 
a'-''-* and b^'^) be 0. Using these parameter values, the initial MA curve is 
assumed to be uncorrelated with duration d and be skewed to the right with 
a peak at a quite early stage (cf.. Figure 2), which is intuitively plausible. 
Second, in (3), we assume that all initial models are linear. Third, to choose 

initial values of (cq^^ , c^°^ ) , for A: = 1 , 2 , . . . , 8 , we first obtain a set of empirical 

estimates oi 7^(^(1 j) at equally spaced time points tj = j x 2/100, for 

A: = 1, 2, . . . , 8 and j = 1, 2, . . . , 49, where tt^'^^p^ (tj^ are defined to be relative 
frequencies of the kth angry behavior in the time intervals [tj — h, tj + h] , and 

h = 0.1 is a bandwidth. Then, (cq^\ c^^^) are chosen to be the estimated least 
squares coefficients of the fitted simple linear regression models from the 
datasets {(MA(tj,a(o),bW),7r^'"'P^(tj)),i = 1,2, . . . ,49}, for /c = 1, 2, . . . , 8. 

By the way, we also tried several different values of Cq'^^^q'''*, and h; algo- 
rithm (15) gives similar results. The iterative algorithm (15) stops at the 

jth iteration when the relative difference between c^"'^ and c^"' ^\ defined 

as ||c^"'^ — c^"' ^ ||/||c^"'^ II , is smaller than or equal to 0.01, where || • || is the 
Euclidean norm. 

To analyze the tantrum data using (l)-(3), we first need to determine 
the orders ma,mb and nik, for k = 1,2, ... ,8, of the polynomial models (2) 
and (3). To this end, the "quasi- likelihood under the independence model 
criterion" (QIC) by Pan (2001) is used. To use this method, we can simply 
treat observations of the eight behavior variables of a given subject at a given 
time point as eight repeated measures of a single binary response. There 
are several slightly different versions of the QIC measure. Here, we use the 
version QIC^ [see Pan (2001) for its definition] that is recommended by Pan 
(2001) and Cui and Qian (2007), due to its simplicity and its high success 
rates in selecting correct models in their numerical studies. To select a final 
model, we start from models (2)-(3) with ma = 2,m6 = 2, and nik = 2, for 
A: = l,2,...,8, which is denoted as "Full" hereafter. Then, a backward model 
selection procedure is implemented, and some related results are summarized 
in Table 2. In the 'Model' columns of the table, we only list ■ma,mf, and m^, 
for = 1, 2, . . . , 8, that are smaller than 2. So, for instance, model '7714 = 
mi = V denotes the one with ma = 2;mb = 2; mi = l;m2 = m^, = 2;m4 = 
1; 771,5 = = "T-y = "T-s = 2. 

From Table 2, it seems that models (2)-(3) with 7774 = 7770 = 777^ = 1 and 
7771 = 7772 = fn^ = 7775 = "^7 = ^-8 = "^b = 2 fit the tantrum data best, by 
the QlCy^ criterion. So, we use these models in the remaining part of data 
analysis. By algorithm (15), the estimated parameter values and standard 
errors (SEs) of these models are presented in Table 3. The SEs are computed 
from the robust variance estimator Vr described in Appendix B. 
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Table 2 

This table lists some models and their QIC ^ values 





Model 




ivioQei 






Full 


12,274.49 








mi — \ 


1 12,304.77 


7714 = Till = 1 


12,349.57 




7712 = '- 


1 12,326.92 


7714 = 7712 = 1 


12,325.00 




7713 = '- 


1 12,302.36 


7714 = 7713 = 1 


12,325.08 




7714 = - 


1 12,272.03 








7715 — '- 


1 12,281.53 


777,4 = 77I5 = 1 






ma = : 


1 12,274.67 


777,4 = 7716 = 1 


12,272.00 




7717 ~ - 


1 12,298.03 


777,4 = 77I7 = 1 


12,301.48 




7718 = : 


1 12,276.45 


7714 = 7718 = 1 


12,275.12 




77i4 = m6=mi = l 12,303.44 


7714 = 7716 = rUa ~ 1 


12,270.78 




7714=7^,6=70,2= 1 12,325.18 


7114 = 7n,6 = 771(, = 1 


12,328.93 




7714 = me = 7713 = 1 12,346.94 


7(14 = 7n,6 = ma ~ mb = 


= 1 12,328.91 








7714 = 7716 = 1, 771a = 


12,330.26 




7714 = 7716 ~ 7fl5 = 1 12,283.75 








7714 = 7716 ~ 7(17 = 1 12,305.91 








7714 = 7716 = 7718 = 1 12,274.97 










Table 3 




Estimated parameters of models (2) and (3) and their standard 


errors (in parentheses) 




when 1714 = 7716 = rUa ~ 1 and mi 


= 7112 = 7713 = 1715 ~ 7717 


— ms ~ mb = 2 


PaTcimeters 


Constcint component 


Linear component 


Quadratic component 


a 




0.0658 (0.0479) 


0.0045 (0.0027) 




b 




0.2410 (0.0725) 


0.0454 (0.0100) 


0.0001 (0.0003) 


Cl 




-6.1432 (0.6562) 


11.3835 (2.5709) 


-8.6876 (2.8639) 


C2 




-1.6569 (0.1318) 


-2.0083 (0.8088) 


3.4437 (1.1259) 


C3 




-1.5895 (0.1338) 


4.0299 (0.7408) 


-4.0841 (1.1209) 


C4 




-5.0790 (0.3770) 


3.1491 (0.6375) 




C6 




-5.2077 (0.4210) 


5.8811 (1.8784) 


-2.6364 (1.9953) 


C6 




-4.5346 (0.3031) 


3.3092 (0.5098) 




cr 




-4.2117 (0.3210) 


6.7421 (1.4772) 


-4.9546 (1.7613) 


Cg 




-3.8145 (0.3129) 


-1.2724 (1.9000) 


2.8053 (2.3235) 



The estimated MA surface is presented in Figure 3(a). From the plot, it 
can be seen that (i) when duration d is smah, MA starts at a quite high 
level and then decreases relatively slowly over time t, (ii) as d increases, MA 
starts at a lower level and drops faster over time, (iii) for a given d value, 
MA peaks at an early time point, and (iv) it seems that the peak of MA 
decreases with duration d when d is small to moderate and then stabilizes 
when d is large. To further demonstrate these results, in Figure 3(b), cross 
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sections of MA when d = 2,5,10,20 and 30 are presented in a single plot, 
from which the above results can be easily seen. It should be mentioned 
that these results are all intuitively reasonable. For instance, in practice, an 
anger episode of a long duration usually starts at a relatively low intensity, 
the peak of its intensity is often at the early stage of the episode, and its 
intensity decreases quite fast in the entire episode. Plots (a) and (b) in 
Figure 3 are made with standardized time t which represents the percentile 
of the duration of a tantrum episode at the corresponding real observation 
time, as discussed in Section 2. For convenience, to perceive the estimated 
MA function in real time, the estimated MA surface and its cross sections 
when d = 2,5,10,20 and 30 are presented again in Figures 3(c) and 3(d), 
respectively, in real time. It seems that all conclusions made from Figures 
3(a) and 3(b) are still true here, except that the real peak time seems to 
increase with duration d when d is small and then stabilizes when d gets 
larger. 

The estimated probabilities of the eight angry behaviors are presented in 
Figure 4 as functions of MA, along with their pointwise 95% confidence in- 
tervals. The confidence intervals are computed based on the robust variance 
estimator Vr. From the figure, it can be seen that (i) Stamp and Throw do 
not happen quite often and their likelihood of occurrence is almost flat over 
the entire range of anger intensity, (ii) Push and Hit do not happen quite 
often either, their likelihood of occurrence increases with MA, and both in- 
creasing trends look close to linear, (iii) the probabilities of Stiffen and Kick 
are at quite low levels in the range of MA, they first increase with MA and 
then stabilize, and (iv) Shout and Scream have relatively large probabilities 
of occurrence, the probability of Shout is quite stable when MA is below 0.4 
and it increases with MA afterward, and the probability of Scream tends to 
increase with MA when MA is below 0.5 and decrease afterward. 

To investigate goodness of fit of the estimated models shown in Figures 3 
and 4, we first apply the Hosmer-Lemeshow test [cf., Hosmer and Lemeshow 
(1989)] to individual angry behaviors. By this approach, observations of the 
kth angry behavior are partitioned into g groups using percentiles of the 
estimated probabilities. Here, we follow the convention that g is chosen to 
be 10 and deciles of the estimated probabilities are used for grouping. Then, 
the Pearson's Chi-square statistic is defined by 

Xi = i:^^V^ for^ = l,2,...,8, 

where 



Ek£= X] 7rfc(tij,a,b,Cfc), 

7rfe(i,j,a,b,Cfe)e{0.1(^-l),0.W] 
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standardized time 



(a) (b) 




real time 

(c) (d) 

Fig. 3. (a) Estimated MA surface is shown in standardized time t. (b) Its cross sections 
when d = 2,5, 10, 20 and 30 are shown by the solid, short-dashed, dotted, dot-dashed and 
long-dashed curves, (c)-(d) Same results as those in plots (a) -(b) are shown in real time. 



7rfc(tij,a,b,Cfe)G(0.1(£-l),0.1£] 

are the expected and observed counts of the ith. group, respectively, and 
ijijk is the observed fcth angry behavior at time Uj. In the above expres- 
sions, when i = l, the grouping interval (0,0.1] can be changed to [0,0.1], 
although this change would not make much difference in the results because 
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Stiffen 



Shout 



Scream 



Stamp 



Push 



Hit 



Kicl< 



Throw 



1 1 1 — 

0.0 0.2 0.4 0.6 



T 1 1 1 — 

0.0 0.2 0.4 0.6 



T 1 1 1 — 

0.0 0.2 0.4 0.6 



T 1 1 1 — 

0.0 0.2 0.4 0.6 



Fig. 4. Estimated probabilities of the eight angry behaviors as functions of MA, along 
with the point-wise 95% confidence intervals for the probabilities. 



Table 4 

Calculated values of X|, for = 1, 2, . . . , 8 



k 


1 


2 


3 


4 




2.03 X 10"" 


2.85 X 10"* 


1.16 


1.63 X 10~^ 


k 


5 


6 


7 


8 




1.53 X 10"^ 


1.38 X 10"^ 


1.59 X 10"^ 


4.76 X 10"^ 



the estimated probabilities would not be exactly in most cases. By these 
formulas, the calculated values of X| are listed in Table 4. 

According to Hosmer and Lemeshow (1989), if the chosen model fits the 
data well, then X| should follow a Chi-square distribution with 8 degrees 
of freedom. The 95th percentile of this distribution is 15.51, which is much 
larger than all values in Table 4. Therefore, the estimated model fits the 
data well in terms of individual angry behaviors. If we would like to study 
how well the estimated models fit the entire dataset, then a reasonable test 
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statistic is 

k=l 

In the two extreme cases that {Xk,k = 1, 2, . . . , 8} are independent of each 
other and that they are perfectly hnearly correlated, the null distribution 
of would be x^(64) and 8x^(8), respectively, with 95th percentiles 83.68 
and 124.06. By Table 4, the calculated value of X'^ from the observed data is 
1.48, which is much smaller than either one of the two 95th percentiles. So, 
by the Hosmer-Lemeshow test, we can conclude that the estimated models 
shown in Figures 3 and 4 fit the observed data well. 

To further check the adequacy of the estimated models, next we use a 
graphical approach by making certain model checking plots, introduced in 
detail in Chapter 22 of Cook and Weisberg (1999). The basic idea of a 
model checking plot is to present the predicted response values based on 
the proposed model and the corresponding empirically estimated response 
values without using the proposed model in a single plot. If the two sets of 
values are close to each other, then we conclude that the proposed model 
describes the observed data well. Otherwise, the proposed model may not be 
appropriate to use. Although this model checking approach is based on our 
visual impression and may not be mathematically rigorous, it is a useful tool, 
especially when related statistical theory is not available. For the tantrum 
data, we can compare the predicted probabilities of the eight angry behaviors 
based on models (l)-(3) and the corresponding empirical estimates of these 
probabilities computed directly from the observed data. Since (i) models 
(l)-(3) assume that these probabilities depend on both duration d and time 
t, (ii) their empirical estimates can be computed at a given time point, 
and (iii) the empirical estimates can not be computed at a given MA level 
due to the fact that MA is unobservable, we choose to present predicted 
probabilities based on models (l)-(3) and the empirical estimates of the 
probabilities when d = 2 or 8 minutes and t = j x 0.05, for j = 1,2, . . . , 19, 
in the model checking plots. The two specific duration values are chosen 
based on the following considerations. First, because most durations in the 
data are quite small (cf.. Figure 1), we can not choose large duration values. 
Otherwise, both empirical probability estimates and predicted probabilities 
based on models (l)-(3) would have large variability. Second, the two chosen 
duration values should be different enough to demonstrate how probabilities 
of angry behaviors depend on duration. The predicted probabilities can be 
computed from the estimated models of (2) and (3). For given values of d and 
t, the empirical probability estimates are taken to be the sample proportions 
of occurrence of the angry behaviors in the time interval [t — 0.05, t + 0.05], 
duration interval [d — 1, d + 1] when d = 2, and duration interval [d — 3, d + 3] 







0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 

time time time time 



Fig. 5. The upper eight panels show the model checking plots of the proposed models 
( 1 )"(3) when d = 2. The lower eight panels show the model checking plots when d = 8. 



when d = 8. The model checking plots for the eight angry behaviors are 
shown in the upper eight panels of Figure 5 when d = 2, and in the lower eight 
panels of Figure 5 when d = 8. Prom the plots, it can be seen that overall 
the predicted probabilities and the empirical estimates of the probabilities 
match reasonably well. For the behavior Scream, the predicted probabilities 
overestimate the trend over time of the empirical probabilities a little bit 
when d = 2 and they match well when d = 8. Small mismatches can also be 
noticed in several other places (e.g., for Hit when d = 2, and for Shout when 
d = 8). 
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4. Summary and concluding remarks. We have presented a procedure 
for modeling both the time course of tantrum anger and the relationship 
between anger and various behaviors that express it. In this model, each of 
eight angry behaviors is assumed to be uniquely driven by a latent variable 
MA that represents momentary anger intensity. MA is modeled by a Beta 
function over the time course of the tantrum. Its two parameters are assumed 
to be polynomial functions of tantrum duration d, and the linkage between 
MA and the angry behaviors is described by generalized polynomial models. 
Model parameters are estimated by a modified version of the GEE iterative 
algorithm, after model selection using the QIC criterion. Goodness of fit of 
the estimated models is checked by the Hosmer-Lemeshow test. Based on 
this analysis, we can make the following conclusions: (i) The selected models 
are appropriate for describing the tantrum data, (ii) The MA function peaks 
near the beginning of the tantrum for a given value of d and then decreases 
gradually. This finding is consistent with the result of "rapid rise and slow 
decline of anger" reported in the psychological literature [cf . , Beck and Fer- 
nandez (1998), Fridja et al. (1991), Tsytsarev and Grodnitzky (1995)]. (iii) 
The peak value of MA appears to decrease with tantrum duration d, in terms 
of the standardized time t. This interesting result was also seen in an earlier 
analysis of a different sort [cf., Potegal, Kosorok and Davidson (1996)]. (iv) 
The probability of the lower angry behaviors Stamp and Throw increases 
only slightly with MA. In contrast, the probability of higher angry behav- 
iors like Shout, Scream, Hit and Kick changes significantly with MA. This 
finding provides an explanation for their differential distribution across the 
tantrum. Namely, the lower angry behaviors are broadly distributed because 
their probability is largely unaffected by changes in anger intensity, while 
the higher anger behaviors tend to appear around the peak of MA. At a 
more fundamental level, this model explains how various behaviors actually 
come to reflect different levels of anger intensity. 

To the best of our knowledge, this analysis is the first to reconstruct the 
time course of anger based on objective behavioral data. It also proposes 
an approach to the analysis of multiple, binary, longitudinal (MBL) data 
with a latent variable involved. Such data are quite common in psycho- 
logical, psychiatric and other areas of social sciences. From a psychological 
perspective, one innovative aspect of the proposed approach is that it starts 
with the temporal relations among various objectively observed behaviors 
in a naturally occurring situation and then works backward to determine 
the corresponding levels of anger intensity, rather than starting with subjec- 
tive estimates of anger intensity and trying to determine properties of the 
associated behaviors. 

At the end of the article, we would like to point out that, besides the 
proposed modeling approach and the corresponding GEE model estimation 
algorithm, there might be other ways to analyze MBL data. One possibility 
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is to use mixed-effects modeling [e.g., Diggle (1988), Laird and Ware (1982), 
Lindstrom and Bates (1988)]. By introducing a random-effects term for de- 
scribing between-subject variability, correlation among repeated measures 
within a subject can be accommodated to a certain degree. It requires much 
future research to set up such a random-effects model in an appropriate 
way, and to compare this alternative approach to the approach discussed in 
this paper. We adopt the GEE approach here because (i) correlation among 
observations does not have to be precisely specified by this approach and es- 
timates are consistent under regularity conditions, and (ii) from a technical 
or mathematical standpoint, it is easier to incorporate the latent variable 
MA in that framework. After the parameters in models (2) and (3) are es- 
timated, it might be of interest to some researchers to compare the effect 
of MA on the eight different angry behaviors. To this end, an appropriate 
multiple comparison procedure should be applied, which also requires much 
future research. 

The tantrum data analyzed in this paper and the R code fitting the final 
model selected by QIC and presented in Figures 3 and 4 and Table 3 are 
available online as supplementary materials [Qiu, Yang and Potegal (2009)]. 

APPENDIX A: PARAMETER IDENTIFIABILITY IN MODELS (l)-(3) 

Theorem A.l. In models (l)-(3), all parameters a,b and c^, fork = 
1, 2, . . . , 8, are identifiable. 

Proof. To verify the identifiability of parameters in models (l)-(3), we 
need to show that, if there are two sets of parameters {sS^\h^^\c^i!'\k = 
1,2,..., 8} and {a'-^\h^'^\cf\k = 1,2, . . . ,8} such that, for all t G [0,1], 

logit(7rfc(t,a«,b«,4^))) 

(5) 

= logit(7rfc(t,a(2),b(2),cf )) for A; =1,2,..., 8, 

then we must have a^^^ = a^^^ , b^-^^ = b^^^ , and c^^^ = c^^-* , for k = 1,2, ... ,8. 

For / = 1,2, let aW = (4'),a«, . . . ,a(il )',bW = (6?\ 6«, . . . , 6« )', 4'^ = 
(Ji) JO JO y 

aW = a« + a?d+ ■■■ + a^lj^^^ , 6« = 6« + h^Pd+ ■■■ + 
Equation (5) implies that, for t E [0, 1], 

c« + cSi) MK{t, a(i) , 6(^)) + ■ • ■ + c« , MA"- {t, a^^) , 6^^)) 

(6) 

= 4? + 4? MA(i, a(^) , 6(^)) + • • • + MA- {t, a^^\h^-)). 
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Obviously, MA(t, a^'), can be written as a polynomial function of t. For a 
given d value, without loss of generality, assume that a^'^ and 6^'^ are positive 
for I = 1,2. Then, the constant terms of the two polynomial functions on the 
two different sides of (6) are c^jj and Cq^-* , respectively. So, we have 

(7) c«=4? forA: = l,2,...,8. 
Equations (6) and (7) imply that, for t G [0,1], 

MA(t, o(i) , ) + ••■ + c^^ MA'"^ {t, o(i) , fo^^)) 

(8) 

= eg) MA{t, a(2), 6(2)) + -.. + c^^l MA™- (t, a('\b^'^). 

The lowest-order terms of the two polynomial functions on the two sides of 

(8) are c^^^t*^" ~^ and cg^t*^" respectively. Therefore, they must be the 
same. Namely, 

cSfc=4fc for A: = 1,2,..., 8 

and 

(9) a(^) = a(2) for all d. 
Using such arguments recursively, we have 



From (9), we have 



c^y=cj^2) for /c = 1,2,..., 8. 



aW=a(2). 



By comparing the highest-order terms of the two polynomial functions on 
the two sides of (8), we have 

(10) = 6(2) for ah d. 

From (10), we have 

b(l)=b(2). 



Therefore, the identifiability of parameters in models (l)-(3) is proved. □ 
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APPENDIX B: MODEL ESTIMATION 



In this part we discuss estimation of the parameters in models (l)-(3) from 
the tantrum data under the framework of generaUzed estimating equations 
(GEE). By this approach, as long as the mean function of observations is 
correctly specified and their variance structure is roughly specified, unknown 
parameters in the models can be estimated by a GEE iterative algorithm 
[cf., Liang and Zeger (1986)]. 

Let yijk denote the binary observation of the kth. angry behavior of the 
ith. child at the time point tij, for j = 1, 2, . . . , rij, i = 1,2, . . . ,N, and k = 
1, 2, . . . , 8, with yjjfc = 1 denoting presence of the behavior and absence. Let 
Yjj = {yiji,yij2, ■ ■ ■ , UijsY be the vector of eight angry behaviors observed at 
the time point tij for the ith child. Its mean vector can be written as 

(11) 7r(%,c) = {iri{tij,a,h,Ci),TT2{tij,a,h,C2),. . . , 7r8(%, a, b, Cg))', 

where c = (a', b', c')' and c = {c[,c'2, . . . , Cg)'. The covariance matrix of Yjj 
can be written as 



where R{tij,a) is a "working" correlation matrix of Yjj which may depend 
on a parameter (vector) a, and 

A{tij,c) = diag(7ri(%,a,b,ci)(l - 7ri(tij,a,b,ci)), 

7r2(%,a,b,C2)(l - 7r2(tjj,a,b,C2)), 

.. .,7r8(%,a,b,C8)(l -7r8(tij,a,b,C8))). 

Then, the generalized estimating equations are defined by 



By Liang and Zeger's (1986) approach, the estimator of c can be computed 
by iterating between the following Fisher scoring algorithm for c and a 

moment estimation for a. Let c'''^^ be an initial estimator of c. Then, in the 
jth iteration, for any j > 0, the updated estimator of c is defined by 



(12) 



V{tij,c, a) = A^/\tij,c)R{Uj,a)A^/\tij,c) 



(13) 




i=ij=i 






(14) 
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Li=ii=i 



X 



dc. 

V (%,c ,a(c ))(Yij -7r(tij,c )) 



where a(c ) is the estimator of a when c is estimated by c 

The "working" correlation matrix can be modeled in several different 
ways. These include the identity matrix Jgxs (when it is reasonable to as- 
sume that the eight angry behavior variables are all independent of each 
other), the compound symmetry structure or, equivalently, the exchange- 
able structure (when it is assumed that correlations between any two angry 
behaviors are all the same), the unstructured pattern where we need to es- 
timate all 8(8 — l)/2 off-diagonal elements of R{tij,a) and so forth. After 
the pattern of R{tij,a) is determined, a can be estimated using the method 
of moments. For instance, if R{tij,a) is assumed to have the compound 
symmetry structure with its {ki,k2)th element being 



Rki ,k2 ) ^) 



1, iiki = k2, 
a, otherwise, 



then a is the correlation coefficient between any two angry behaviors and 
its moment estimator is given by 

^ i=lj=lki<k2 

where A^* = ^fLi Sj4i(8 x 7)/2, p is the number of parameters in the mean 

function [cf., equation (11)], njk = {yijk - '^{tij,c^ ^^))/[7v{tij,c''^ ^^)(1 - 

Tv{tij,c'^ ''))]^^^5 s-iid c'^ ^ is the estimated parameter vector obtained 
from the GEE algorithm [cf., equation (14)]. Obviously, in the above ex- 

pression, {rj^fc} are the standardized residuals and a(c ) is the pooled 
sample correlation constructed from these residuals. 

Let c be the solution of (13). According to Liang and Zeger (1986), under 
some mild regularity conditions, A^/^(c — c) is asymptotically Normal with 
mean zero and covariance matrix 

^«^.L"^|:|(^)V-'fe...)(^?^))^' 



d7v{tij,c) 
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X Cov(Yij)T/^H%,c,a) 



d7v{tij,c) 

dc 




,c,a 




where Cov(Yjj) denotes the covariance matrix of Yjj. In apphcations, Vr 
can be estimated by its finite-sample version after c and a are replaced 
respectively by c and a, and Cov(Y.y) by (Y^- — 7r(tjj, c))(Yjj — 7r(tjj,c))'. 
The resulting estimator Vr is often called the "robust" variance estimator. 

A favorable property of the GEE method is that the above-mentioned 
asymptotic normality of c depends only on the correct specification of the 
mean function [cf., (11)]; it does not depend on the correct choice of the 
"working" correlation matrix R{tij,a) in (12). For this reason, if we do 
not have any prior information about the correlation matrix of Yjj, then 
R{tij,a) is often chosen to be the identity matrix in applications [cf., Diggle 
et al. (1994), Section 8.4.2]. 

In models (11) and (12), besides a, there are 30 parameters in the pa- 
rameter vector c. In the updating formula (14) of the conventional GEE 
algorithm, we need to compute the inverse of the 30 x 30 matrix 



N rii 

VV(a7r(t,,,e^^"'V5c)V-i(t 



i=ij=i 



)/dc) 



in each iteration, which requires a great amount of CPU time. To reduce the 
computing burden, we suggest dividing the vector c into g blocks, that is, 
c = (C]^, 62, . . . , c^)', and then using the following Fisher scoring algorithm 
with blocked parameters. In the jth iteration, for any j > 0, update the 
estimator of , for ^ = 1, 2, . . . , g, successively by the updating formula 

-0-1) 



(15) 



+ 



N rii 

EE 

U=lj=l 



d7v{tij,C 



dcf 



■ N Tli 

EE 

i=ii=i 



F^H^.e*, «(£*)) 



97r(t 



dCe 



d7T{tij,C ) 



where c 



(ci ,. 



V ^{tij,c,a{c)){Yij - 7r{tij,c)) 
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SUPPLEMENTARY MATERIAL 

Tantrum data and R Code (DOL 10.1214/09- AOAS242SUPP; .zip). This 
is the tantrum anger data analyzed in the paper. The data has 10 columns. 
The first 8 columns denote the binary status of the 8 angry behaviors, with 
1 denoting "present" and "absent." The 9th column is the duration of a 
tantrum episode, and the 10th column is the standardized observation time. 
The data are ordered by duration (i.e., the 9th column). This is a R code 
fitting the final model selected by QIC presented in Figures 3 and 4 and 
Table 3 of the paper. 
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